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Long data outages may occur in applications of global navigation satellite 
system technology to orbit determination for missions that spend significant 
fractions of their orbits above the navigation satellite const ellation(s). Cur- 
rent clock error models based on the random walk idealization may not be 
suitable in these circumstances, since the covariance of the clock errors may 
become large enough to overflow flight computer arithmetic. A model that is 
stable, but which approximates the existing models over short time horizons 
is desirable. A coupled first- and second-order Gauss-Markov process is such 
a model. 

INTRODUCTION 

According to Van Dierendonck, 1 the Allan variance technique represents the power spec- 
tral density, £(/), of a clock’s frequency fluctuation as a function of frequency by 

S(f) = h 2 f 2 + hi/ + ho + h-i/- 1 + h_ 2 /- 2 (1) 

where h^ is the “white phase noise,” hi is the “flicker phase noise,” ho is the “white frequency 
noise,” /i_i is the “flicker frequency noise,” and h- 2 is the “random walk frequency noise.” 
The random walk model developed in Brown and Hwang 2 makes use of the frequency 
noise terms only, i.e. /io, h- 1, and h- 2. These terms give rise to asymptotes on an Allan 
variance plot, an example of which Figure 1 shows. The white frequency noise asymptote is 
yJho/2r ) and dominates at small averaging times. The flicker frequency noise asymptote is 
>/(2 ln2)/i_i, and dominates at intermediate averaging times. The random walk frequency 
noise asymptote is 27 ry/^ 2^r/6, and dominates at larger averaging times. Figure 1 also 
shows the Allan variance of the random walk model developed in Ref. 2, and the Allan 
variance of a Gauss-Markov model, discussed below (the appendix of this work specifies the 
random walk model) . 

As Figure 1 shows, the random walk model nicely approximates the asymptotic Allan 
variance plot over all time scales, while avoiding certain model representation issues related 
to the flicker frequency noise. 2 The possible difficulty with this model is evident in Figure 2, 
i.e. the clock bias error envelope becomes larger than the diameter of the earth within two 
days. Nominally, global navigation satellite system receivers process measurements from 
four or more satellites at intervals on the order of seconds, so the error growth Figure 2 
shows does not occur. However, in a medium, highly-elliptical, or geosynchronous Earth 
orbit mission, even nominal satellite visibility conditions may not permit rectification of the 
time bias estimate for long periods of time. Allowing the uncertainty to continue to grow 
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Allan Variance 



Figure 1 The Gauss-Markov model closely approximates the Allan 
variance of the random walk model over a wide range of sample times. 



Time History 



Figure 2 The random walk model is unstable; the Gauss-Markov model is stable. 
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as Figure 2 shows appears harmless enough, since a Kalman filter should be able to quickly 
re-converge to reasonable uncertainty levels once visibility conditions permit additional 
measurement processing. In fact, a very large covariance is somewhat desirable in such 
cases, if the filter has a residual edit function that could be triggered by large residuals that 
might occur at measurement re- acquisition. 

As stated at the outset, the desirability of this approach will be negated if the covariance 
becomes large enough to overflow the flight computer’s arithmetic, causing the filter software 
to crash. Figure 2 indicates that, for the somewhat average clock parameters simulated here, 
this condition probably would not occur for several days. (This statement assumes the use 
a U — D factorized filter, which stores quantities related to the square root of the covariance; 
this problem would be accentuated in a filter storing the total covariance matrix.) Since 
the random walk model is much simpler than the Gauss-Markov model, in most cases it 
may be adequate. 

If on the other hand a requirement exists for a filter to operate through very long outages, 
or if a filter is tuned or mechanized in such a way as to create the possibility of covariance 
overflow, a coupled first- and second-order Gauss-Markov model proposed herein may be 
used. As Figures 1-3 show, the model is stable, reaching a specified steady-state value after 
a few hours, and it closely approximates both the time histories and the Allan variances of 
the random walk model over shorter time intervals. Note that we do not mean to imply 
that real clocks have bounded error growth; we merely suggest that open-loop-stable error 
models may be more suitable for navigation filters in some applications. 




Figure 3 Over short intervals, the Gauss-Markov and random walk 
models agree well. 
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PROBLEM STATEMENT AND SOLUTION 


t 


The Gauss-Markov clock bias model proposed herein is a coupling of a first-order Gauss- 
Markov (FOGM) process driving the bias state directly, and a second-order Gauss-Markov 
(SOGM) process driving the bias state via integration of its direct effect on the clock drift. 
Let b represent the clock bias, and d represent the clock drift. Then, the continuous-time 
representation of the model is 


b(t) 


— 1/r 1 

' b(t) ‘ 

I 

Wl(t) ‘ 

d(t) 


-“>n -2C Un . 

_ d(t) _ 

I 

. w 2 {t) _ 


where r is the FOGM time constant, w\ is the FOGM zero-mean driving white noise 
with power spectral density (PSD) <71, uo n is the SOGM natural frequency, Q is the SOGM 
damping ratio, and is the SOGM zero-mean driving white noise with PSD q\ 2. Note that 
the SOGM is not merely an integrated FOGM; as the appendix shows, the latter is not 
stable. 


Although Eq. 2 is a stochastic differential equation, we can seek a mean-square solution. 
With x = [b, d }' and w = [w\, w 2 ]\ Eq. 2 becomes: 


x{t) = Ax(t ) + w(t) 


( 3 ) 


where A is given by 


A = 


— 1/r 1 

~ u} n ~‘^C UJ n 


( 4 ) 


Since E [w] = 0, the mean solution, #(£), will be given by the homogeneous solution of Eq. 2. 
The covariance, P = E [(# — x) (x — ^) 7 ] , where 


P = 

Pll 

V 12 



Pbb 

Pbd 


_ P21 

P22 


_ Pbd 

Pdd _ 


must satisfy 

P(t) = AP(t) + P(t)A' + Q, 

where Q = diag(<?i, q 2 ). 


( 5 ) 

( 6 ) 


Analytical Solution 

Let us first find the homogeneous solution for Eq. 2, in the form: 

x(t) = $(t,to)x(to). (7) 

The state transition function #(t, t G ) must possess the following properties: 

#(Mo) = A$(t,t 0 ), $(t 0 ,t 0 ) = I 

= $(t,ti)$(ti,t 0 ), t 0 < t\ < t 
t Q ) = <f>(A t), At = t - t 0 
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( 8 ) 

( 9 ) 

( 10 ) 



A state transition function that satisfies these properties is the following: 
p aAt 

*(A t) = — 

where 


p = 1/ T > 

(12) 

a = — 2 (P + 2£a> n ) , 

(13) 

b = + PC^n ~ -^P 2 1 

(14) 

Ud = - c 2 > 

(15) 


b cos bAt + (a + 2 (,uj n ) sin bAt sin bAt 

b cos bAt + (a + (3) sin bAt 


-<n 2 sin bAt 


( 11 ) 


and we assume that b 2 > 0. Since the noise input is zero-mean, this state transition matrix 
allows one to find the mean solution, x(t), to Eq. 2. 

Let 

P,, 

c — — 2" + C 

then, the covariance is given by the following: 

„2 


pn(At) = qi 
+ 


f e 2aAt _ i 


4a 

e 2aAt cos: 


c“ \ g 2 aAt g j n 2b At ( b 2 — c 2 + 2ac 

^ b 2 / 4(a 2 + b 2 ) 


2bAt — 1 / ab 2 — ac 2 + 2b 2 c 


4(a 2 + b 2 ) 


b 2 


/e 2aAt - 
+ b 2 V 4a 


1 e 2aAt sin 2bAt + a cos 2b At) — a 


p 2 2 (At) = 92 


a 2aAt 


- 1 / <r \ e 

( 1+ p' + 


2aAt 


4 (a 2 + b 2 ) 
sin 2bAt ( b 2 — c 2 + 2ac 


4(a 2 + b 2 ) 
e 2 aAt cos — If ab 2 — ac 2 + 2b 2 c 


4(a 2 + b 2 ) 

9i< 4 /e 2aA t _ ! 


+ 


b 2 


b 2 


e 2aAt (b sin 2bAt + a cos 2b At) — a 


4a 


(16) 


(17) 


Pi 2 (At) = 


b 2 


+ b 2 


f (i-e 

4a v 


2a _|_ 


A (!_ e 2aAt ) + 
4 a v 7 


4(a 2 + b 2 ) 

p 2aAt ^p c _ a ^ g j n 2&At + (ac — b 2 ) cos 2bAt] — (ac — b 2 ) 

4(a 2 + b 2 ) 

p 2aAt _|_ fc c ) s j n 2bAt + (ac — b 2 ) cos 2bAi] — (ac — b 2 ) 


4(a 2 + b 2 ) 


(18) 


Approximate Solution 

In online filtering applications, the analytical solution may prove cumbersome to imple- 
ment, and additionally, the time step may be small in comparison to the modes of the 
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k 


The formal solution to 

(19) 

( 20 ) 

DISCUSSION 

Characteristics of the Solution 

Examining the solution given above, we see that the parameter a governs the rate of 
decay of all of the exponential terms. Therefore, we define the “rise time” as that interval 
within which the transient response of the covariance will reach a close approximation to 
the above steady-state value; thus, we define the rise time as follows: 

tr = (21) 

a 

Next, we note that all of the trigonometric terms are modulated by 2 b; thus we may view 
this value as a characteristic damped frequency of the coupled system. The period of the 
oscillation, II, is then 

n = t r/b (22) 

In the limit as t — > oo, all the exponential terms in the analytical solution die out, so that 
the steady-state value of the covariance simplifies to: 

p/„\ = 1 [ <?2 + (2a 2 + b 2 + c 2 - 2ac)qi fi(q 2 + <?iu> 2 ) 

v ~ 4a(a 2 + b 2 ) [ (S(q 2 + <?iw 2 ) (2a 2 + b 2 + c 2 + 2 ac)q 2 + qiuj* 

(23) 


system. In such cases, an approximate solution may prove useful. 
Eq. 2 is 

P{t + At) = 3>(Af)P(f)$'(Af) + S(At), 


where 


4>(r)Q4> / (r)dr, 

the simplest approximation to which is S' (At) = QAt. 


Comparison of Numerical, Analytical, and Approximate Solutions 

Figure 4 compares the approximate numerical and the exact analytical solutions for the 
time evolution of the covariances, for a particular set of paremeters. Initially, the position 
error is a few percent, but it quickly decreases to a tiny fraction of a percent. The velocity 
error is never more than a fraction of a percent. In this example, the period of the oscillatory 
response, II, is about 8.7 hours. 


Parameter Sensitivity 

Figure 5 shows the results of a sensitivity analysis. The left and right columns show the 
sensitivity of the formal error in bias and drift, respectively, to parameter variations. The 
FOGM time constant varies from half a day to over a year, the SOGM damping varies from 
a very underdamped value of 10 -5 to critical, and we varied the SOGM natural frequency 
so as to cover a similar range of bias and drift values to the ranges the other parameter 
variations produced. It is notable that varying uj n produces a somewhat wider range of 
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Figure 4 The approximate analytical solution agrees well with the exact solution. 


variability in bias than it does in drift, relative to the ranges of variation that the other 
parameters produce. This is similar to the relationship between bias and drift error that 
would result from a purely SOGM model, <Jdl&b — which the appendix shows. For this 
reason, we find that uj n is a useful parameter for adjusting the relative magnitudes of bias 
and drift in steady-state, as the discussion below describes further. 


Tuning 

To use the coupled FOGM/SOGM process as a clock model, we wish to tune it such 
that it approximates the usual RW model over time intervals associated with filtering. For 
a FOGM, choosing the time constant, r = 1 //?, to be much larger than the sampling 
interval causes the FOGM to resemble a RW over intervals less than the time constant. 
A similar effect occurs for a SOGM, by choosing the product £uj n much smaller than the 
sample interval. The corresponding parameter for the coupled FOGM/SOGM process is 
1/a (recall that a = —/3/2 — (u n ); choosing 1/a to be much larger than the sampling interval 
will cause the coupled FOGM/SOGM process to resemble a RW over intervals shorter than 
1/a. It may also be desirable to select a value of £ near unity so as to minimize any ringing 
that otherwise might creep into the coupled model from the SOGM dynamics; however, we 
have not found this to be a significant problem for even very underdamped values, as the 
sensitivity analysis described above shows. 

Next, we wish for the slopes of the covariances of the coupled model to be similar to the 
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Bias Sensitivity to x Drift Sensitivity to x 
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Drift Sensitivity to £ 



Figure 5 Sensitivity of formal error to parameter variations around 
baseline of r = 86400 sec, u n = 0.0001 rad/sec, £ = 0.075009, 
qi = 0.017 m 2 /sec, q 2 = 0.027 m 2 /sec 3 . 


8 







RW model. The slopes for the coupled model are given by 


P22(t) = <?2 - 2ujlp 12 (t) - 4Cw„P22(i) (24) 

Pi2(t) = P22(t) ~ (P + 2Qu n )p\ 2 (t) - ulpn(t) (25) 

pn(t) = 2pi 2 (t) + q\- 2f3pi\(t) (26) 

whereas for the RW clock model, the slopes are given by 

P22 = q2 (27) 

Pi 2 (t) = p 22 (t ) = q 2 t (28) 

puit) = 2pi 2 (t) + qi = q 2 t 2 + qi (29) 


Given that /3 and £ u n have been chosen to be small*, we can make the slopes approximately 
equal over time intervals less than 1 /a by choosing q\ and q 2 to be the same as in the RW 
model. 

The final step in tuning the coupled model is to check the steady-state values to ensure the 
magnitudes are acceptable. We have found that adjustment of u n is an effective means for 
modulating the relative size of the steady-state drift versus the steady-state bias variance. 


Ensembles of Sample Realizations 

Figures 6 and 7 show a couple of examples of ensembles of sample realizations for two dif- 
ferent parameter selections. These examples use the same process noise PSDs, but different 
values of the FOGM/SOGM parameters. Both use rather underdamped values of C- The 
product is of the same order in both examples however. The FOGM time constant is 
almost an order of magnitude larger, and u n is an order of magnitude larger, in Fig. 6. In 
both cases, the rise time is about 30 hours. The most notable difference between the two 
examples is the behavior of the sample realizations of the bias error, with the bias ensemble 
in Fig. 7 appearing less noise-like, which may or may not be of interest to filter performance. 

CONCLUSION 

The coupled first- and second-order Gauss-Markov model for clock errors which this work 
contributes may prove useful in a variety of applications. The model is open-loop stable, 
which means the designer of a navigation filter may definitively specify a maximum value for 
the magnitude of the error growth in the absence of measurements, thereby alleviating non- 
sensical results and numerical overflow issues. The model closely approximates a random 
walk over a designated rise time, so that over a given time interval of interest, the model 
may be tuned to match the error characteristics of real clocks whose random errors are 
adequately modeled by the Allan variance characteristics described in Refs. 1 and 2. If the 
analytical solution is not suitable for implementation in some cases, a remarkably simple 
approximation for the process noise covariance contribution over small time steps has shown 
to be adequate. 

* Since C is on the order of 1, then terms multiplied by ufa may also be neglected. 
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Figure 6 Sample Realization 1. 
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Figure 7 Sample Realization 2. 
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APPENDIX 

Here, we specify the random walk, integrated first-order Gauss-Markov, and second-order 
Gauss-Markov models. 


Random Walk Model 

The random walk clock model of Ref. 2 is given by 


■ bit) ' 


' 01 ' 

' fe(t) ' 

+ 

Wl(t) 

. - 


0 0 

. d(f) . 

_ w 2 (t) _ 


which leads to the following state transition matrix, 


$(A t) = 


1 At 
0 1 


and process noise covariance, 


S(At) = 


q 1 At + q 2 ^f q2^f 

n A{2 

Q2 — 


q 2 At 


(30) 


(31) 


(32) 


Clearly, the clock drift variance increases linearly with elapsed time, and the clock bias 
increases as the cube of elapsed time. 


Integrated First-Order Gauss-Markov Model 

The integrated first-order Gauss-Markov model is given by 


b(t) 


' 0 

1 

' b(t ) 

+ 

0 

. . 


_ 0 

-!/r _ 

. d (*) . 

. w 2 (t) _ 


which leads to the following state transition matrix, 


#(Af) = 


1 r( 1 — e- At / r ) 
0 g-At/i- 




(33) 


(34) 
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and process noise covariance, 

~ r 2 { (! _ e -2A t/r)2 + m + 4 (1 _ e-At/r) | r ( X _ g-At/ry 


5(At) = 


<?2T 


(1 - e~ At / T ) 2 (1 - e~ 2At / T ) 


(35) 


Clearly, this is an unstable model, as the clock bias variance increases linearly with elapsed 
time. 


Second-Order Gauss-Markov Process 

The second-order Gauss-Markov model is given by 


1 

1 

- 

1 

l 

L 


0 1 
~ u) n ~^‘C UJ n 


b(t) 

d{t) 


+ 


0 

m{t) 


Reference 3 gives the following solution for the second-order Gauss-Markov process. 

g ^UJnt g C^nt 

E \b(t) ] =b 0 (co d cos co d t + (u n sin u>dt) + d 0 sinu> d t 

Wd w d 

e -(Unt e ~CUnt 


E \d(t) ] = — b 0 uJri sin Udt + d Q (ud cose — C Un sina;^) 

™d ™d 

g 2£c o n t 

1 j — (^d + ^C^n^d cos u d t sin ujfjt + 2£ 2 u; 2 sin 2 u d t) 


E [& 2 (t)] 


E [d\t)} = 


02 

4C ^ 
02 
4C 


1 - 


w d 
-2< u>„t 


- (u; d - 2C ) uj n uj d cos LU d t sin u d t + 2( 2 u; 2 sin 2 ui d t) 


Wa 


E [b(t)d(t)] =^2 e 2< ’ W " i sin2 Udt 


(36) 

(37) 

(38) 

(39) 

(40) 

(41) 


In the over-damped case, replace sin and cos with sinh and cosh, respectively. In the 
critically-damped case, let u d — > 0, ( —> ► 1: 


E [6(i)] =b 0 e Unt (l + a > n t) + d Q te Wnt 
E [ d(t )] = — b 0 uj 2 te~ u ’ nt + d 0 e~“ nt ( 1 - u > n t) 

E i b2 (*)] =£3 [! " e _2Wnt (l + 2u; n t + 2u; 2 t 2 )] 


E [d 2 (t)] =-^- [1 - e- 2 ^'(l - 2w„t + 2 u> 2 n t 2 )\ 

E[6(t)d(t)] =^ e - 2bJnt 


In any case, as t — > 00 , 


P(t — > 00 ) 




4C w n 


iA4 0 

0 1 


Thus, the ratio of the steady-state standard deviations of x and x will be 

°d _ 

— ^71) 

^6 


(42) 

(43) 

(44) 

(45) 

(46) 

(47) 

(48) 


12 



and these are related to the power spectral density by 


